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O ABSTRACT 

\Q ■ 

o : 

Motivated by the presence of numerous dark matter clumps in the Milky 
Way's halo as expected from the cold dark matter cosmological model, we conduct 
Oh 1 numerical simulations to examine the heating of the disk. We construct an initial 

galaxy model in equilibrium, with a stable thin disk. The disk interacts with dark 
matter clumps for about 5 Gyr. Three physical effects are examined : first the 
mass spectrum of the dark matter clumps, second the initial thickness of the 
galactic disk, and third the spatial distribution of the clumps. We find that the 
^ ■ massive end of the mass spectrum determines the amount of disk heating. Thicker 

disks suffer less heating. There is a certain thickness at which the heating owing to 
the interaction with the clumps becomes saturates. We also find that the heating 
produced by the model which mimics the distribution found in Standard CDM 
cosmology is significant and too high to explain the observational constraints. 
On the other hand, our model that corresponds to the clump distribution in a 
ACDM cosmology produces no significant heating. This result suggests that the 
ACDM cosmology is preferable with respect to the Standard CDM cosmology in 
explaining the thickness of the Milky Way. 

Subject headings: cosmology — dark matter — galaxies:kinematics and dynamics 
— methods: numerical 
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1. Introduction 

Hierarchical clustering governed by cold dark matter (CDM) is widely believed as a 
cosmological scenario which is responsible for the growth of the structures in the universe. 
According to the hierarchical scenario, small dark matter halos should collapse earlier, but 
later fall into larger structures. The process of smaller halos being assembled into a larger 
halo does not always destroy the smaller ones, thus hierarchical structures are seen in many 
objects, such as clusters of galaxies. 

Recent high-resolution simulations have successfully shown that hundreds of galaxy-size 
DM halos survive in clusters of galaxies (Ghigna et al. 1998; Colin et al. 1999; Klypin et al. 
1999a). A remarkable outcome of the high-resolution cosmological simulation in Standard 
CDM model by Moore et al. (1999) even shows that survival of substructures or satellites 
occurs not only on cluster scales, but also on galactic scales. They show that a galaxy 
of a similar size as the Milky Way should contain about 500 satellites, which is, however, 
much more than the number of the observationally identified satellites. That many satellites 
should survive, was confirmed also by Klypin et al. (1999b) and Ghigna et al. (2000). Klypin 
et al. note that the results of the Standard CDM simulation are close to those of a ACDM 
simulation with the same circular velocity function of substructures. This indicates that the 
prediction of the presence of many satellites is a general outcome of the hierarchical scenario 
and not particularly dependent on the cosmological models. 

Compared with the observational results, these cosmological models yield a large number 
of the satellites, approximately a factor of 50 more than the number of satellites observed in 
the vicinity of the Milky Way. 

Klypin et al. (1999b) suggest that the problem of the missing satellites could be resolved 
(i) by identification of some satellites with the high-velocity clouds observed in the Local 
Group (Blitz et al. 1999) or (ii) by considering dark satellites that failed to gather enough 
gas to form stars, because of expulsion of gas the supernova-driven winds or because of gas 
heating by the intergalactic ionizing background. The latter possibility implies that the halos 
of galaxies may contain substantial substructures in form of numerous invisible DM clumps. 

A statistic of strong gravitational lensing is an approach to identify the dark clumps 
in the Milky Way. Chiba (2002) indeed finds evidence for the existence of the numerous 
satellites in the Milky Way. 

If a great amount of the dark satellites exist within the Milky Way's halo, their inter- 
action with the disk might cause disk heating. In their Standard CDM model, Moore et al. 
(1999) found that a large fraction of the substructures have very eccentric orbits, so that 
they could cause resonant heating of the disk, and even heating by impulsive shocking owing 
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to their penetration through the disk. 

On the other hand, it is known that the Milky Way has a quite thin disk, whose scale 
height is about 200 pc. From the "thinness" Toth and Ostriker (1992) have made an ener- 
getical analysis of the disk heating owing to accretion of matter, and derived the constrain 
that the Milky Way should have accumulated no more than 10 % of the disk mass within 
the past 5 Gyr. This estimation of the disk heating might, however, be too large because 
the actual interaction between the disk and satellites is more complicated. For example, 
a single satellite could dissolve before reaching the disk (Huang and Carlberg 1997) or the 
energy injected into the disk could excitate coherent warping motions rather than heat the 
disk (Velazquez and White 1999). For the case of the interaction of the disk with many 
substructures in the halo, additional detailed numerical investigations are necessary. Font 
et al. (2001) have studied the case of the ACDM cosmological model, and found that the 
substructures are not efficient perturbers for heating of the disk, since the masses of the 
clumps are lower than those of the clumps predicted in the Standard CDM model, and also 
because the clumps are located far away from the disk and seldom get near the disk. 

The disk kinematics and dynamics is a good probe not only for examining the cosmo- 
logical models, but much more generally for clarifying the halo substructure that is difficult 
to derive from direct observations. Therefore in this study we aim to investigate disk heating 
by dark matter (DM) clumps for a wide range of parameters. 

Numerical experiments on the disk dynamics are not an easy task, especially when 
studying the vertical structure, because of the wide range of the dynamical scales among 
the different components in galaxies. The smallest scale is disk, with a scale height about 
200 pc, while the largest scale is the dark halo extent of > 100 kpc. Many numerical studies 
consider disks of 700 pc (0.2 x the disk scale length) in thickness (Velazquez and White 1999; 
Font et al. 2001). The question is whether the heating rate obtained for such thick disks is 
the same as for thin disks like the real Galactic disk. 

In this study we construct our initial galactic models following Kuijken and Dubinski 
(1995), which is nearly in exact equilibrium and which allows us to set up disks as thin as 
the real Galactic disk. Several additional observational constraints are taking into account 
to build the galactic models. We study three physical effects that could affect disk heating; 
first the mass spectrum of the DM clumps, second the initial thickness of the galactic disks, 
and the third the spatial distribution of the clumps. 

This paper is organized as follows. In section 2 we present the galaxy model which 
provides a very stable thin disk comparable to the real Milky Way disk. A model of a 
clumpy dark matter halo is presented in section 3. Numerical simulations of the interaction 
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between disk and the clumps are specified in section 4. Section 5 presents the results of our 
numerical simulations on the disk heating by examining the effects on the mass spectrum 
of the clumps, the initial disk thickness, and the spatial distribution of the clumps. We 
summarize and discuss our results in section 6. 



We use the self-consistent disk-bulge-halo galaxy model given by Kuijken and Dubinski 
(1995, hereafter KD). This model provides nearly exact solutions of the collisionless Boltz- 
mann and the Poisson equations, so that one can control subtle differences in the initial 
conditions. This is important especially to set up very thin equilibrium disks, which is the 
most difficult part. Using this model we could successfully realize an equilibrium disk with 
realistic thickness of only 200 pc. 

Here we summarize the model properties briefly. The configurations of the disk, bulge 
and halo are determined by their distribution functions. The bulge distribution function 
depends only on the relative energy E, while that of the halo depends also on the angular 
momentum about the symmetric axis, L z . 

The halo distribution function takes the form 



The relative energy E is defined so that E — at the edge of the distribution, where 
p = 0. The halo distribution function has five free parameters: the potential at the center 
^/q, which appears implicitly in the definition of the relative energy, the velocity scale (Jo, 
and three factors A, B, and C . The factors A and B control the system flattening (q) and 
the core radius (R c ), respectively, and all the three factors are scaled by the density scale 



The bulge distribution function takes the same form of a King model (Binney and 
Tremaine 198 7) which is given by 



2. Galaxy Model 




[(AL 2 Z + B) exp(-E/a 2 ) + C] [exp(-E/a 2 ) - 1] HE < 0, 
otherwise, 



(1) 



(Pi)- 




p b (27r^)-3/2exp[(* - * C )K] {exp[-(£ - * c )/o-g] - 1} 



if E < * c , 

otherwise. 



(2) 



This distribution depends on three parameters: the cutoff potential \& c , the center density 
Pb, and the velocity dispersion a h . 
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For halo and bulge distributions, the densities are given by analytic functions of R and 

*. 

The disk distribution function must depend on three integrals of motion so that it can 
sustain a triaxial velocity ellipsoid within the disk. Since the third integral is not analytic, 
KD employ an approximate form of the distribution function which depends on E, L z and 
the vertical energy E z = ^(R, z) — ^(R, 0) + \v 2 z . E z is only an approximate integral, but 
well conserved for stars in nearly circular orbits in a warm disk. This quantity is the only 
approximation in KD's model. The distribution function takes the form 



f ( F T F \ — Q ( Rc) P ~ dW r,-n 

/diskl P ' " z) (2n^(R c ) a R \R c )~a z {R c ) 6XP 



E p — E C (R C ) E z 



(3) 



where E p = E — E z , R c and E c are the radius and energy of a circular orbit with angular 
momentum L z . Q and k are the circular and epicyclic frequencies at radius R c , respectively. 
The 'tilde' functions pj, <7r, and a z are free functions. Therefore the density and the radial 
velocity dispersion are conveniently selected as 
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* 2 (fl,3z d ) 



(4) 



and 

4 = 4(0)exp(-i?/ J R d ). (5) 

Pd and a z are iteratively adjusted so that the densities on the mid-plane and at height z = 3z d 
will agree with those given by eq. (4). The distribution of the disk has six free parameters: 
the disk mass Md, the radial scale length i?d, the vertical scale height z&, the disk truncation 
radius -R ou t) the truncation width <5-R out , and the central velocity dispersion of the disk cxr(0). 

The choice of the above parameters is made to satisfy observational properties of the 
Milky Way (see also Tsuchiya (2002)), that is the solar radius Ro = 8 kpc, the circular 
velocity of the disk at the solar radius V c = 220 km/s, the total surface density within 1.1 
kpc of the disk plane S L1 (_R ) = 71±6M pc~ 2 (Kuijken and Gilmore 1991), the contribution 
of the disk material to S L1 , Ed(i?o) = 48 ± 9M Q pc~ 2 (Kuijken and Gilmore 1991), the total 
Galaxy mass within 50 kpc M tot (< 50kpc) = 5.4l°; 2 x 10 n M o (Wilkinson and Evans 1999), 
and the total Galaxy mass within 170 kpc M tot (< 170kpc) = 1.9l?;f x 1O 12 M (Wilkinson 
and Evans 1999). 

We assume that the disk mass is M d = 5 x 1O 1O M , the disk scale length R d = 3.5 kpc, 
the disk scale height z& = 245 pc, and the disk truncation radius and the truncation width 
are R ont = 7i?d and 5R = i?d, respectively. The bulge mass is about 15 % of that of the disk, 
while the tidal radius is 2.38 kpc. The contributions of the bulge and disk to the rotation 
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curve is shown in Fig. 1. We also examined thicker disks with = 525 pc and 1050 pc, as 
described in the next section. 

For the halo, we assume A=0 in eq. 1, which means that the halo is nearly spherical, 
except for the flatness caused by the presence of the disk potential. As for the extent of the 
halo, we consider two different models. One has the virial radius at 262 kpc, with a dark 
matter mass of 8.59 x lO n M . We refer to it as the standard halo model. The other has 
the virial radius at 1.35 Mpc, with a mass of 4.1 x 10 12 M Q , which is called the extended 
halo model. The difference between the two models lies only in the outer part of the halos. 
The inner density profiles (R < 50 kpc) are about the same (Fig. 2 and Fig. 4), thus the 
rotation curves of the bulge, disk, and the halo are nearly equivalent between two models 
(Fig. 1 and Fig. 3). The extended halo model is introduced only to produce different clump 
distributions as described below. 

The halo and the bulge are treated as an external force, while the distribution of the 
disk is realized by self-consistent particles. 

Without the substructure clumps in the halo, this model is fairly stable. The net increase 
in the disk thickness is only due to 2-body relaxation among the disk particles, which is only 
Azd ~ 150 pc after 4.76 Gyr. 



3. Dark Matter Clumps Model 

The numerical cosmological simulations (Moore et al. 1999) have demonstrated the 
presence of the dark matter clumps within dark halos. We distribute the clumps in a similar 
way as shown in the numerical results. 

The dark matter clumps are represented by rigid bodies with a Navarro, Frenk & White's 
density profile (Navarro et al. 1997). It provides an accurate fit to the density profiles of 
CDM halos: 

Pcrit^c 

(r/r B )(l + r/r s y 

Here p crit = 3H 2 /8nG is the critical density for closure of the universe, and 5 C and r s are the 
characteristic density and the scale radius of the clump. In practice, it is more convenient 
to use mass and concentration to characterize each clump instead of 5 C and r s . 



P(r) = ; ; ° ; , 2 - (6) 



The mass of a dark matter clump, M, is related to its virial radius r 20 o, which is defined 
as the radius within which the average density equals 200 times the critical density at the 
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present p C rit,o, by the relation 



M 



200 



200p crit ,o- 



(7) 



3 



Then the characteristic density 5 C is expressed in terms of the concentration c = r 2 oo/r s as 



Each clump has two characteristic parameters, but we fix the concentration to a typical 
value c = 20. Since the interaction between the clumps and the disk is dominated by distant 
encounters, the adopted value for the concentration does not affect the heating process so 
much. 

The simulations of Klypin et al. (1999b) and Ghigna et al. (2000) show that the total 
mass of the clumps is about one tenth of the halo mass. We therefore give the clumps 10 % 
of the halo mass of the standard model, which is 8.59 x 1O 1O M in total. The distribution 
of positions and velocities of the clumps is assumed to be the same as the halo distribution 
function (eq.l), so that the clump distribution does not change in time. This is the most 
reasonable assumption which is important in order to maintain an equilibrium in the whole 
system. Otherwise the change in the mean potential of the halo might cause undesired 
change in the disk thickness. The background halo component, which is treated as an fixed 
external field, is reduced by 10 %. 

The orbital properties of the clumps could be represented by the radii of their pericenters 
and apocenters. Fig. 5 shows the distribution of the pericenters and the apocenters of the 
clumps. The lines of the constant eccentricity, which is defined as e = (r a — r p ) / (r a + r p ) 
are superposed in the plot. The eccentricity distribution has a median at 0.6, which is in a 
good agreement with the Standard CDM simulations that predicts high eccentricity orbits 
with a median apocentric-to-pericentric distance of 6:1 (Moore et.al, 1999). 

For the mass spectrum of the clumps, we examine three different possibilities. Model 1 
and model 3 have only one single population of clumps with clump masses of 10 8 M Q and 
1O 9 M , respectively. These two models are meant to be extreme cases. Model 2 has a mass 
spectrum with a mass range of 10 8 to 1O 9 M according to the probability distribution 



following the cosmological prediction by Moore et al. (1999), Klypin et al. (1999b) and 
Ghigna et al. (2000). This model is considered to be the most realistic examination. In this 
model the masses are assigned randomly to the clumps. These three models were designed 
in order to examine the role of the individual mass of the clumps on the disk heating. 



200 



(8) 



3 ln(l + c) -c/(l + c)' 



N{M)dM oc M~ 2 dM, 



(9) 
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Two additional physical effects on disk heating are studied. One is the initial disk 
thickness. Thinner disks have smaller vertical velocity dispersions, they should therefore be 
more sensitive to disturbances from the clumps, and might experience a large thickening. 
Therefore we ran models 4 an d 5 which have an initial disk scale height of 525 pc and 
1050 pc, respectively. Both models have the same mass spectrum as model 2. 

Another influential property of the clumps on disk heating is the position of their peri- 
center, where the gravitational tides become the largest. The probability distribution of the 
pericentric radii of the DM clumps is shown in Fig. 6. The abscissa shows the pericentric 
radii of the clumps, and the ordinate shows the cumulative fraction of the number of clumps 
which have pericenters less than the value on the abscissa. About 18 % - 20 % of the clumps 
have pericenters within the solar radius. This is similar to the clump distribution obtained 
by Moore et al. (1999) in the Standard CDM cosmology. 

In order to examine the effect of the position distribution of the clumps, we created 
model 6, which has the same mass spectrum as model 2, but the Milky Way has the density 
profile, according to the extended halo model. This model has two noticeable features; (i) 
the inner density profile is the same as the standard halo model so that all the kinematical 
properties of the disk are equivalent, and (ii) the tidal radius is 5 times larger than that 
in the standard model. We distribute the clumps in the same way as the halo particles. 
Thus their spatial distribution is also 5 times larger in radius than in the other models. In 
this model, the total clump mass is the same as before. Therefore, the mass fraction of the 
clumps to the whole halo is 2 %. The probability distribution of the pericenters is shown 
in Fig. 6. The number of clumps that have pericenters within the solar radius is about 3.5 
%. This model corresponds to the clump distribution predicted from a ACDM cosmology, 
as demonstrated by Font et al. (2001). 

The characteristic parameters of all models are summarized in Table 1. 

4. Numerical Simulation 

We calculate the interactions between the disk particles by using a hierarchical tree 
algorithm Hernquist (1987) with a tolerance parameter of 9 to \ = 0.7. The time integration is 
made with the leap-frog method and a fixed time step of 1.75 Myr. The softening length of 
the disk particles is 70 pc. We use 131,072 particles in the disk. The forces from the clumps 
are calculated by direct summation. We have followed the evolution up to 4.76 Gyr. 

The calculations are made on Pentium III - Linux workstations (800MHz). The typical 
calculation time is 90 seconds per time step. 
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5. Results 

In this section we show the results of our numerical simulations on the disk heating. 
Three different physical effects are examined here: the mass spectrum of the clumps, the 
initial disk thickness, and the spatial distribution of the clumps. 

Throughout the section, we are measuring the disk heating by the increase of the disk 
thickness, Az d , and by the change in the velocity dispersion, A<jr, Aa z . To calculate these 
quantities, the disk is stratified in concentric annuli, each of which containing 8192 particles. 
The disk thickness is evaluated simply by the mean square deviation of the ^-coordinate of 
disk particles from the center of the disk plane, z d (R) = {z d 2 ) l l 2 , within the annuli. 

Even without the clumps, the self-heating within the disk owing to two-body relaxation 
among the disk particles is always taking place in numerical simulations with modest number 
of particles. In our simulations with 131,072 particles in the disk, the increase owing to the 
self-heating is typically Az d = 150 pc, and Aa z = 4 km/s at the solar radius after 4.76 Gyr 
simulations. These values are reasonably small, but need to be taken into account when 
analyzing the results. The curves corresponding to the self-heating after 4.76 Gyr are shown 
as internal heating in Fig. 7 - Fig. 13 . 

5.1. Mass Spectrum of the Clumps 

The results of the disk heating obtained from the simulations of the model 1, 2, and 3 
are shown in Fig. 7 and Fig. 8. The comparison between the results makes the effect of the 
mass spectrum of the DM clumps clearer. 

The growth of the disk thickness in model 1 is Az d ~ 615 pc, 630 pc, and 900 pc at 
R — 2, 8, and 16 kpc, respectively, while those values in model 3 are 860 pc, 1298 pc, and 
2197 pc. It is clear that the interaction of a galactic disk with a few, but massive clumps 
(model 3) is more effective to heat the disk than the interaction with many but less massive 
clumps. This result can be explained by the estimate in an impulsive approximation, which 
predicts that the increase in the energy of the disk is AE oc N M 2 . If there is a mass 
spectrum, dN(M) oc M~ 2 dM = M _1 d(logM), AE is dominated by massive clumps. The 
disk thickening obtained from model 2 (Fig. 7b) is, in fact, much closer to that in model 3. 
This result could be also understood by the impulsive approximation, as the energy input 
to the disk from the logarithmic interval in the mass spectrum is dominated by the more 
massive clumps; dAE oc MdilogM). Therefore the higher limit of the mass spectrum has 
an essential contribution on the disk heating. 
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A difference feature between models 2 and 3 is their dependence of the disk heating on 
the disk radius. For all models the outer part of the disk is heated up more than the inner 
part, but this gradient is the largest in the model 2. The increase of the disk thickness for 
R < 13 kpc is smaller in model 2 than in model 3, while it is opposite in the outer disk. In 
model 2, which imitates the results of the Standard CDM model, the disk thickness becomes 
1 kpc at the solar radius after 4.76 Gyr interactions. 

An alternative way to estimate the disk heating is the increase in the velocity dispersion. 
Fig. 8 shows the same results as Fig. 7. The existence of a few massive clumps has a large 
influence on the disk heating. This figure shows the increase in the radial (the left panels) 
and the vertical (the right panels) velocity dispersions. After 4.76 Gyr, the interaction of 
the disk with clumps in model 2 produces a R = 60 km/s and a z = 35 km/s. The properties 
of disk heating in model 2 seem to agree with those of the thick disk component at the 
solar radius which has <jr = 63 km/s and a z = 38 km/s and a vertical scale length of 1 kpc 
(Freeman 1996). 

It is noticeable that the initial velocity dispersion is not isotropic or ctr is larger compared 
to a z . After 4.76 Gyr, in the inner part of the disk (R < 5 kpc), there is no increase in 
the radial velocity dispersion while the vertical velocity dispersion increases a lot. This 
anisotropy in the heating is seen also in the outer part of the disk in model 2. A possible 
explanation is that the heating owing to the clumps is aiming at a more isotropic velocity 
dispersion. 

5.2. Initial Disk Thickness 

Our model 2 has an initial disk thickness comparable to the observed value in the Milky 
Way, which is about 250 pc. In many simulations in other works, however, the typical initial 
disk thickness was chosen much larger. For example, Font et al. (2001) employ the initial 
disk with a thickness of 700 pc. The susceptibility of the disks to the heating owing to the 
interaction with the DM clumps should depend on the initial disk thickness. Thicker disks 
are possibly less sensitive and suffer less heating. In order to clarify the difference in the 
disk heating according to the initial disk thickness, we carried out three investigations with 
the models 2, 4, and 5. These models have almost the same bulge, halo, and clump systems, 
but the initial disk thicknesses are 245 pc, 525 pc, and 1050 pc, respectively. The results of 
the disk heating are shown in Fig. 9, and Fig. 11. 

After 4.76 Gyr, the disk in model 2 achieves a thickness of 1 kpc at the solar radius, 
thus the increase rate is about 400%. In models 4 and 6, the final disk thicknesses at the 
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solar radius are 1.1 kpc and 1.5 kpc, respectively. The increase rates in the thickness are 
thus 210% and 70%. This tendency is more pronounced in the inner regions. In the central 
region R < 3 kpc, the disk in model 5 has experienced nearly no heating. On the other hand, 
in the outer regions R > 14 kpc, the heating owing to the interaction with the clumps is so 
large that the difference in the initial disk thickness is negligible. This feature is more clearly 
seen in Fig. 10. The panels in the figure show the time variations of the disk thickness (a) at 
R = 2 kpc, (b) at the solar radii, and (c) at R = 16 kpc, respectively. It is natural that the 
outer part of the disk experiences more heating, because it is closer to the DM clumps. A 
common characteristic to all radii is that the difference in the initial disk thickness decreases 
in time. The curves in the panel (a) indicate that the heating will be saturated at the level 
of Zd ~ 1.2 kpc. In the outer parts the heating is still going on to make the disk even thicker. 
The growth of the velocity dispersions in the disks, which is shown in Fig. 11, also gives the 
same information namely that thicker disks suffer less heating. 

This flaring effect, that is the increase in the disk thickness especially in the outer part, 
is seen in VLA observations of HI disk of galaxies NGC 100 and UGC 9242 (Bosma 1991). 
The interaction with the hypothetical DM clumps would be an explanation of this feature. 

5.3. Spatial Distribution of the Clumps 

In our models the clump's spatial distribution is the same as the halo distribution in 
order to ensure the equilibrium. This clump distribution, however, creates more clumps in the 
inner part of the halo than that obtained by the cosmological simulations. For example, Font 
et al (2001) employ a clump distribution in which only 3 % of the clumps have pericenters 
within the solar radius. Our model 6 is designed to remove the clumps from the inner 
region. This is accomplished by introducing a more extended equilibrium halo model. We 
distribute the same clumps as in the model 2 in the extended halo. The spatial extent of 
the clump distribution in this model becomes 5 times larger, and only 3.5 % of the clumps 
have pericenter within the solar radius (see Fig. 6). 

The increase in disk thickness and velocity dispersion is shown in Fig. 12 and Fig. 13. 
In the model 6, the disk heating owing to the interaction with the clumps is comparable 
with the disk internal heating, up to the disk edge. In this case, we would need more 
precise simulations with larger number of particles, but it is already clear that even with the 
presence of the many clumps surviving in the halo, tidal disk heating is negligible and the 
disk remains as thin as the real Milky Way disk for about 5 Gyr. This result is in agreement 
with conclusion of Font et al. (2001). 
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6. Discussion and Conclusions 

We have made a series of numerical simulations, examining heating of galactic disks 
that are embedded in a halo with many small DM clumps, as suggested by cosmological 
simulations. We have built up the initial galaxy model as precisely in equilibrium as possible, 
so that we could simulate a stable disk that is as thin as the Milky Way. 

We have shown that with the mass spectrum N(M)dM oc M -<I.\I. that is consistent 
with cosmological simulations, the massive end of the mass spectrum determines the amount 
of disk heating. The number of the clumps at the massive end could be only a few, so that 
the massive limit of the spectrum might fluctuate among galactic halos. This fact suggests 
a variety of disk thickness because of the random fluctuation in the number of the most 
massive clumps. 

As a second result, the simulations demonstrate that thicker disks suffer less heating. 
It is possible that a disk with 700 pc thickness shows no significant heating whereas another 
disk with 200 pc thickness suffers significant heating in the same circumstances. This also 
means that there is a critical thickness at which the heating owing to the interaction with 
the clumps is saturated. 

Finally, we have considered the relation to cosmology. Our model 2 mimics the clump 
distribution found in the Standard CDM cosmology (Moore et al. 1999; Klypin et al. 1999b). 
The disk heating by the clumps has already been discussed by Moore et al. (1999). They 
argued by a simple impulsive approximation that disk heating would be too high to explain 
the age-temperature relation for disk stars (Wielen 1974). Our numerical simulations are 
supplementary to their analysis, and we confirm that the heating is significantly high. The 
thickness of 1 kpc at the solar radii is comparable to the present thick disk, but this heating 
takes place within 5 Gyr and even a young thin disk would have no chance to survive. On 
the other hand, our model 6 corresponds to the clump distribution of the ACDM cosmology, 
similar to the one Font et al. (2001) have studied, but more closely representing the Milky 
Way. As found by Font et al. (2001), this clump distribution causes no significant heating 
of the disk. These results suggest that the ACDM cosmology is preferable with respect to 
the Standard CDM cosmology in explaining the thickness of the Milky Way. 

We thank Burkhard Fuchs and Andreas Just for their useful comments and suggestions. 
EA is grateful to Astronomisches Rechen-Institut (ARI), Heidelberg, in particular Rainer 
Spurzem, for allowing her to conduct simulations on ARI's Unix cluster. TT is indebted to 
Alexander von Humboldt foundation for research grant. 
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Fig. 1. — Circular velocity profiles for the Galaxy model with the standard halo. The cross 
sign (+) shows the observational circular velocity at the solar radius. 
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Fig. 3. — Circular velocity profiles for the extended halo model. The cross sign (+) shows 
the observational circular velocity at the solar radius. 
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Fig. 4. — Volume density profiles of the components in the extended halo model. 
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Fig. 5. — Distribution of the pericenter and the apocenter of the clumps in model 1. 
cross (+) symbols indicate the individual clumps. The lines of constant eccentricity e = 
0.4, 0.6, and 0.8 are also shown. 
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Fig. 6. — Cumulative distribution of the pericenters of the clumps in the models 1, 2, 3, and 
6. The ordinate shows the fraction of the clumps which have a pericenter smaller than the 
value on the abscissa. In the models 1 to 3, 20 % of the clumps have pericenters smaller 
than 10 kpc, while the corresponding fraction is only 5 % in model 6. 
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Fig. 7. — The growth of the disk vertical scale length after interaction with DM clumps with 
different masses given in (a) model 1 (b) model 2 (c) model 3. 
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Fig. 8. — Disk velocity dispersion in vertical and radial direction after 4.76 Gyr of interaction 
with different mass distributions of DM clumps given in (a) model 1 (b) model 2 (c) model 3. 
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Fig. 9. — The growth of the disk vertical scale length which initially has height (a) 245 pc 
(b) 525 pc (c) 1050 pc, after interaction with DM clumps, specified in model 2, model 4 and 
model 5. 
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Fig. 10. — Heating of a disk with initial height 245 pc, 525 pc and 1050 pc (a) at it! = 2 kpc, 
(b) at the solar radius, and (c) at R = 16 kpc. 
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Fig. 11. — Velocity dispersion of a disk which initially has a height of (a) 245 pc (b) 525 pc 
(c) 1050 pc, after 4.76 Gyr of interaction with DM clumps, specified in model 2, model 4 
and model 5. 
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Fig. 12. — The growth of the disk vertical scale length produced by (a) model 2 where 18 % 
of the clumps have pericenter radii less than the solar radius and (b) model 6 with only 3.5 
% of the clumps crossing the disk within the solar radius. 
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Fig. 13. — Disk velocity dispersion profiles produced by models with a fraction of clump 
pericenter radii inside the solar radii of (a) 18 % in model 2 and (b) 3.5 % in model 6. 
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Table 1. DM clumps and disk models. 



Model 


Mass of each clump 


Number of clumps 


Disk height 


Model 1 


10 s M 


859 


245 pc 


Model 2 


10 8 - 10 9 M 


335 


245 pc 


Model 3 


10 9 M Q 


86 


245 pc 


Model 4 


10 8 - 10 9 Mq 


335 


525 pc 


Model 5 


10 8 - 10 9 Mq 


335 


1050 pc 


Model 6 


10 8 - 10 9 Mq 


315 


245 pc 



